Influence of correlations on the velocity statistics of scalar granular gases 
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The free evolution of inelastic particles in one dimension is studied by means of Molecular 
Dynamics (MD), of an inelastic pseudo-Maxwell model and of a lattice model, with emphasis on the 
role of spatial correlations. We present an exact solution of the Id granular pseudo-Maxwell model 
for the scaling distribution of velocities and discuss how this model fails to describe correctly the 
homogeneous cooling stage of the Id granular gas. Embedding the pseudo-Maxwell gas on a lattice 
(hence allowing for the onset of spatial correlations), we find a much better agreement with the MD 
simulations even in the inhomogeneous regime. This is seen by comparing the velocity distributions, 
the velocity profiles and the structure factors of the velocity field. 

I. INTRODUCTION 

Some twenty years ago, S.Ulam proposed an elegant and simple model Q aimed to show how the Maxwell distri- 
bution for the velocity statistics in a gas may be derived directly from the conservation laws encoded in the particle 
dynamics. He considered the evolution of an assembly of N particles, whose initial velocity distribution is arbitrary. 
At each instant a random pair of particles is selected and their velocities are updated as if they had collided elastically. 
Such a process converges spontaneously to a macro-state whose velocity distribution is a Gaussian, with the variance 
given by the constant energy of the system. 

Relaxing the energy conservation law, i.e. considering inelastic collisions between particles, leads to different 
statistical properties, which in fact are related to a well defined physical problem, i.e. the dynamics of Fluidized 
Granular Materials g). These are assemblies of moving macroscopic particles, whose evolution is controlled by the 
inelasticity of their mutual collisions. They present unusual structural and statistical properties which stimulated many 
studies not only because of their industrial and technological relevance, but also due to the fundamental questions 
they raise. 

Roughly speaking, the inelasticity causes two phenomena on the velocity statistics: a) in the absence of energy 
injection the fluid cools down since the variance of the velocity distribution, the so called granular temperature, 
decreases; b) the form of the distribution may even cease to be Gaussian and approach a stationary shape under a 
suitable rescaling. 

A large amount of recent studies || indicates that the evolution of a freely evolving granular system consists of three 
stages: a homogeneous cooling during which both the velocity and the density are spatially uniform; a second stage 
where the velocity field forms vortices and a third stage where the density becomes highly inhomogeneous and clusters 
appear. Unfortunately, many of the existing approaches have intrinsic limitations: MD numerical simulations are time 
demanding and hardly access to the correlated regimes; kinetic theories, based on the molecular chaos assumption, 
seem to be limited to the quasi elastic situation; hydrodynamic equations suffer from their phenomenological deriva- 
tion. For these reason minimalistic models, which either lend themselves to analytical solutions or greatly reduce the 
computational effort, can provide useful hints. 

Hereafter we shall focus on the physics of rapid granular flows employing a one dimensional collision rule in order 
to render the analysis as simple as possible. This choice, in spite of its apparent simplicity shares many features 
with the physics of higher dimensions. The paper is organised as follows: first, we shall revisit the inelastic hard rod 
model moving on a ring by means of event driven dynamics. Next, we shall compare these results with those of the 
inelastic Ulam model (which corresponds to a gas of pseudo-Maxwell molecules Jll],[l3|]), for which we find the exact 
asymptotic velocity statistics, showing how they disagree even in the homogeneous regime. Thirdly, we shall extend 
Ulam's model including spatial fluctuations, and show that many features of the real one-dimensional gas, both in 
the homogeneous and in the inhomogeneous regime, are recovered. 



II. INELASTIC HARD RODS ON A RING 



We define a one dimensional granular gas as a system of hard rods, confined to a ring and colliding inelastically. 
After a binary collision the scalar velocities of the particles change according to: 
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where r is a restitution coefficient, which takes on the value f for perfectly elastic systems and for completely 
inelastic particles. Few years ago McNamara and Young and Sela and Goldhirsh ||, simulated such a model 
and observed a universal algebraic decay of the kinetic energy, E(t) oc t~ 2 (Haff's Law), together with an anomalous 
behavior for the global velocity distribution, even in the early homogeneous regime. Furthermore, the appearance of 
strong inhomogeneities is the precursor of a numerically catastrophic event, named inelastic collapse: i.e. particles 
perform an infinite number of collisions in finite time interval. A renewal of interest on one-dimensional granular flows, 
has been generated by the recent work of Ben-Nairn et al. 0. They eliminated the inconvenience of the inelastic 
collapse, by means of a physically motivated 'regularization': they assumed the restitution coefficient to be elastic for 
small momentum transfer. This allows the system to enter into a new dynamical regime, independent of the choice 
of the regularization. During such a regime, the kinetic energy decays as E(t) oc i~ 2 / 3 , for any r < 1. A direct 
inspection of the hydrodinamic profiles, shows that such a regime is highly inhomogeneous, with density clusters and 
shocks in the velocity field. They suggest that inelastic systems behave asymptotically as a sticky (r = 0) gas M, 
which is known to be described by the Burgers equation in the inviscid limit g. This should justify the presence of 
shocks and clustering and, moreover, predicts that the tails of the velocity distribution decay as exp(—t; 3 ) ||. On the 
other hand, an accurate numerical test of this prediction j(| seems problematic, because the bulk of the distribution 
does not deviate appreciably from a Gaussian (see Fig. ^). Such a behavior reflects the fact that the asymptotics is 
dominated by the dynamics of clusters of particles, which move through the system and coalesce, like sticky objects. 
Although the sticky nature of the inelastic granular dynamics seems a very appealing description, it is not clear if it 
applies to higher dimensions p7[. 
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FIG. 1. Rescaled velocity distributions for the MD and in the lattice gas, during the homogeneous. The initial distribution 
(both models) is also shown. The distributions refer to systems having the same energy. Data refer to N = 10 6 (both models) 
particles with r = 0.99 and r — 0.5 (for the lattice model in the inhomogeneous regime). 



To the best of our knowledge, exact analytic treatments of the hard rod model are limited to the homogeneous 
regime, by means of solutions of the Boltzmann equation. Caglioti et al. [[[o| solved asymptotically the Boltzmann 
equation for a Id inelastic hard rod gas in the limit r — ► l~ and p — > oo (where p is the density). They obtained 
the asymptotic velocity distribution, which, if rescaled to have a constant variance, is the superposition of two delta 
functions (see als ojflq ] ) . For MD simulations, the observation of two distinct peaks in the rescaled velocity distribution 
was reported in [gp(see Fig. |l|): it represents a precursor of that singular asymptotic distribution. However, the 
appearance of inhomogeneities and, eventually, the inelastic collapse prevent the approach to the limiting distribution. 
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FIG. 2. Rescaled velocity distributions for the ID MD and in the ID lattice gas, during the inhomogeneous phase. The 
distributions refer to systems having the same energy. Data refer to N = 10 6 (both models) particles with r = 0.99 and r = 0.5 
(for the lattice model in the inhomogeneous regime). 




FIG. 3. Asymptotic velocity distributions P(v,t) versus v/vo(t) for different values of r from the simulation of the inelastic 
pseudo-Maxwell (Ulam's) model. The asymptotic distribution is independent of r and collapse to the Eq.(^). The chosen initial 
distribution is drawn (same result with uniform and Gaussian initial distribution), more than N = 10 particles. 



III. ULAM'S MODEL 



Alternatively, one may attack the Boltzmann equation, by assuming a simpler form for the scattering cross section 
in the collision integral, i.e. taking the relative velocity of the colliding pair to be proportional to the average thermal 
velocity: \v — v'\ ~ \[E, This is the so called pseudo-Maxwell model |13|| . The energy factor \f~E can be eliminated 
via a time reparametrization, and one obtains a simpler equation: 

d T P(v, t)+P(v, t)=P jdu P{u, t)P (/3v+(1-/3)u, t) (3) 

where = 2/(l+r) and the t counts the number of collisions per particle. Eq. (||) is the master equation of the inelastic 
version of Ulam's scalar model: at each step an arbitrary pair is selected and the scalar velocities are transformed 
according to the rule of Eq. (||). In a stimulating paper, Ben-Nairn and Krapivsky considered such scalar model, 
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and obtained the evolution of the moments of the velocity distributions. Since at large times, (v n ) ~ exp(— ra n ), 
and the decay rates a n ^ na,2/2 (they depend non-linearly on n), they argued that such a multiscaling behavior 
prevents the existence of a rescaled asymptotic distribution / such that P(v,t) — > f {v / vq(t)) / vq{t) , for large r, 
where Vq(t) — J v 2 P(v,T)dv — E(t). On the contrary, we believe that such "multiscaling" behavior only indicates 
the fact that the moments of the rescaled distribution J x n f{x)dx — (v n ) /vq diverge asymptotically for n > 3, and 
does not rule out the possibility of the existence of an asymptotic distribution with power law tails. In fact, the 
Fourier transform of Eq. (|[) 
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FIG. 4. Structure factors S(k, t) against kt 2 ^ 3 for the ID MD and against fcr 1//2 for the ID lattice gas model, in the 
inhomogeneous phase. Fimes are chosen so that the two systems have the same energies. Data refers to system with more than 
N = 10 5 particles, r = 0.5 (both models). 
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possesses several self-similar solutions of the kind P(k,r) = /(fcwo(r)), which correspond to the asymptotic rescaled 
distribution P(v,t) = f (v / vq(t)) / vq(t) . Many of them do not correspond to physically acceptable velocity dis- 
tributions p| . The divergence of the higher moments implies a non analytic structure of / in k = 0, since 
(v n ) /vq — (— i) n -§rnf{k)\h=Qi and represents a guide in the selection of the physical solution. As shown in Fig. [| our 
data collapse on the function 
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corresponding to the self-similar solution f(k) = (1 + |fc|) exp(— \k\) Notice that (H) is a solution of Eq.(0) for 

every r < 1, i.e. the asymptotic velocity distribution does not depend on the value of r < 1, as shown in Fig. pT(left 
frame). However, such solution (||) does not resemble the distribution obtained in the simulation of the inelastic hard 
rod model, suggesting that the approach d la Ulam fails to describe even the homogeneous cooling phase (see Fig. |l|). 
We believe that such failure is due to the absence of velocity correlations. 



IV. ONE-DIMENSIONAL LATTICE MODEL 



To reinstate spatial correlations we considered a Id lattice model where the N sites have a "velocity" Vi. At every 
step of the evolution a pair of neighboring sites is chosen randomly and undergoes an inelastic collision according to 
Eq.(||) if (vi — vj) x (i — j) < 0. The latter condition, which avoids collisions between particles "moving" away from 
each other, is to be be referred below as the kinematic constraint and plays a key role in the formation of structures 
during the inhomogeneous phase. A unit of time r correspond to N of such steps (i.e. to 1 collision per particle on 
average). Note that, since there is no real particle movement, the number of collisions per particle, r, remains the 
only measure of time in our model. This represents a problem in the direct comparison with MD, because r(t) is a 
well defined "universal" function for the homogeneous regime only: in the inhomogeneous regime r(i) diverges due to 
the inelastic collapse, or depends on the "regularization" of the restitution coefficient. We observe that initially the 
energy E(t) decreases exponentially, as in the Haff regime for MD simulation, but after the formation of long range 
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velocity correlations, the decay slows down, and a power law decay of the form E ~ r -1 / 2 , appears. As stated above, 
this law cannot be directly compared to the t~ 2 / 3 behavior, observed in MD, since, in the non-homogeneous regime 
the function r(t) depends on the regularization employed to avoid the inelastic collapse. However the problem can be 
circumvented if one assumes the energy, instead of r, as the physical clock of both dynamics. 

A surprising agreement between the MD and the lattice gas is found in the comparison of the velocity distribution 
taken at instants corresponding to the same energy. Both the early two humps structure and the late Gaussian 
shape seen in MD simulations agree with those of the lattice gas, as shown in Figs. [I] and g. The Gaussianity of the 
asymptotic field distribution, and the r -1 / 2 asymptotic decay suggest that, in the lattice model, the scalar velocity 
field v(i, t) is governed by a discretized version of the diffusion equation of the form: d T v(i, r) = vAv(i, r), where A 
is a lattice laplacian. We checked this possibility studying the structure factor S(k, t) = v(k, r)v(—k, r) where i(k, r) 
is the Fourier transform of the field v(i,r). In fact, at late times, S(k,r) = fiikr ' 2 ) scales similarly to a diffusive 
process as shown by the data collapse in Fig. || (right frame). However the form of the scaling function differs from 
the Gaussian shape predicted by the diffusion equation. In same figure (left) we report the structure factor of the MD 
simulation (fourier transform of C(r, t) — (v(i, t)v(i + r, t)), where v(ijt) is the velocity of the i-th particle, see below). 
In this case a good collapse, S(k,t) = f2(kt 2 / 3 ), has been obtained J20|. The two collapses can be directly compared, 
because, in both models, the structure factor is a function of k/E, consistently with the previous discussion on t and 
r. The algebraic tails of the structure factors observed in both models carry important information about the nature 
of the growing velocity correlations. A comparison of velocity profiles at instants having the same energy in the two 
models is presented in Fig. ||. The top frame displays v(x,t) of the MD, with the characteristic shock and preshock 
structures. In the middle frame, it is shown the profile v(i, t) = v(xi(t)) of the MD, where i is the particle label. The 
difference between v(x, t) and v(i, t) in MD is in the sign of large gradients: v(x, t) displays large negative gradients, 
while v(i, t) has large positive gradients. Finally, the bottom frame displays the v(i, r) profile of the lattice gas which 
fairly compares with the analogous MD quantity. The power law behavior S(k) ~ k~ 2 observed in both models are 
due to the presence of short wavelength defects, viz. shocks, as predicted by Porod's law [|l9j. Note that the presence 
of shocks in the lattice model is due to the kinematic constraint, since they disappear, together with the Porod's tails, 
relaxing the constraint. 



10000 

0.04 



-0.04 
0.04 



11000 



particle position 

12000 13000 14000 



-0.04 
10000 



11000 



11000 



12000 13000 

particle number 

12000 13000 

I I 1 1 



15000 



-J 1 1 

: shock p 


1 i 

•e-shock 









2 

re 
n 
c 



3 

a 



14000 



14000 




15000 
15000 



o 

Q. 



FIG. 5. Portions of the instantaneous velocity profiles for the ID MD (top v(x,t), middle v(i,t)) and for the ID lattice gas 
model (bottom, v(i,r)). In the middle frame we display the MD profile against the particle label in order to compare the 
shocks and preshocks structures with the lattice gas model (the dotted lines show how shocks and preshocks transform in the 
two representations for the MD). Data refers to TV = 2 ■ 10 4 particles, r — 0.99 (both models). 



V. CONCLUSIONS 



In summary, we have investigated the velocity statistics of the homogeneous and inhomogeneous cooling regimes 
of the Id granular gas. The complex observed behavior, featuring two distinct dynamical regimes, the onset of 
long-ranged velocity correlations and shocks in the velocity profiles can be reproduced even disregarding the real 
kinematics of the particles. In fact we introduced a lattice model of immobile particles which reproduces quantitavely 
the velocity statistics both at short and long times. The structure factors and the velocity profiles can also be 
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qualitatively compared. The crucial ingredient of the model is the onset of velocity correlations in the system. Such 
correlations are relevant for the velocity statistics even when short-ranged. We corroborated such a conclusion through 
the study of the mean-field limit of the lattice model, which corresponds to a pseudo-Maxwell model of an inelastic 
gas. In this limit the velocity distributions change drammatically, displaying large tails, at odds with the real gas. 
We found the exact expression for this power law tailed distribution, which represent a self-similar solution of a non- 
trivial integro-differential Boltzmann like equation |lrj| . Although the one-dimensional character of our results may 
seems restrictive, they can be interesting in several contests, as kinetic theory and microscopic models for turbulence 
(Burgers equation) or for surface growth (KPZ equation). 
We wish to thank E. Caglioti and A. Vulpiani for many fruitful discussions. 
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